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Abstract 

An accurate method to compute the annihilation rate in positronic systems by means of 
quantum Monte Carlo simulations is tested and compared with previously proposed methods 
using simple model systems. This method can be applied within all the quantum Monte Carlo 
techniques, just requiring to accumulate the positron-electron distribution function. The annihi- 
lation rate of e + LiH as a function of the internuclear distance is studied using a model potential 
approach to eliminate the core electrons of Li, and explicitly correlated wave functions to deal 
with all the remaining particles. These results allow us to compute vibrationally averaged an- 
nihilation rates, and to understand the effect of the Li + electric field on positron and electron 
distributions. 

PACS number(s): 36.10.-k, 02.70.Lq 
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I. INTRODUCTION 



In positron and positronium (Ps) chemistry and physics, the annihilation rate T 27 plays an 
important role since it correlates with many aspects of the local environment where the positron 
annihilates. For instance, "pick-off" annihilation in solutions and in solid materials, "on the fly" 
annihilation in atomic and molecular gases, and bound state annihilation of positronic compounds 
are just few of the experiments where T2 7 can be measured and successively interpreted. 

Although these experiments are relevant both technologically and scientifically jl|, ||, only few 
theoretical studies have been devoted to accurately compute annihilation rates for realistic systems 
like atoms and molecules in order to compare with experimental data or to predict trends (jjj, |], ^ 
® @' H' [| 0i 0' Moreover, these have been restricted to deal at most with four active electrons, 
so that only a bunch of systems have been studied so far. We believe this scarceness of results to be 
primarily due to the intrinsic difficulty in obtaining accurate wave functions for larger systems, and 
to the computational effort requested with respect to ordinary matter compounds when standard 
ab initio methods are employed . 

For these reasons, quantum Monte Carlo (QMC) methods Jll| represent an alluring alternative 
to these methods, to Density Functional Theory, and to explicitly correlated wave functions in 
computing energies and properties of realistic positronic systems. QMC techniques are well described 
in the literature, so we avoid to burden this paper with the details of the methods and constrain 
ourselves to discuss only the technical issues relevant for the specific problem. 

Not requiring the analytical calculation of integrals, QMC allows one to use any physically 
sensible wave function. This possibility increases the chances to obtain an accurate description of any 
class of systems once all the relevant physical information is included in the chosen analytical form of 
the wave function. Having defined a trail wave function *S?t for a system, QMC techniques allow one 
to compute the differential and non-differential properties of the system by sampling ^t^o, or 
vPq. Here, is the exact ground state function of the system. This task is usually accomplished by 
creating a distribution of points (also known as configurations or walkers) in configurational space 
whose density is proportional to the aforementioned ^>^, ^fx^o, or 'J'q. 

Keeping in mind the above remarks, it might appear that the QMC methods should accurately 
predict any interesting observable for positronic systems. This is indeed correct except for extremely 
local operators like Dirac's delta 8, and hence for T2 7 that is proportional to its expectation value, 
for which an accurate sampling of small configurational space volumes is needed. These operators 
are well known to represent a challenge for QMC due to the discrete nature of the configuration 
ensemble and the finite length of the simulations. 

As far as the mean value of the Dirac's delta 8 operator is concerned, one faces an additional 
difficulty when trying to estimate its mean value. Even admitting a perfect sampling in the regions 
where two particles are close to each other, the primitive method of counting the number of times 
the interparticle distance r is smaller than a given radius r w (i.e. counting the ones that fall into 
a spherical well of radius r w ) || has an associated statistical error that diverges for r w — > 
This fact means that the estimation of the statistical error of the extrapolated value is based on 
shaky grounds. 

Although not a solution, a slightly better approach was devised by substituting the simple sphere 
with a Gaussian function centered at the coalescence point H, |[ [Tl[] . The variance of this estimator 
also goes to infinity upon decreasing of the Gaussian width, but it diverges less faster than the 
one of the spherical well, therefore allowing a statistically more accurate estimation of (8(r [_)) = 

Due to the interest in computing (5(r)} for many systems, attention has been paid to solve these 
problems, and remedies have been suggested in the framework of all the QMC techniques. 

As far as variational Monte Carlo (VMC) is concerned, different methods have been proposed 
that may solve this difficulty allowing one to compute the needed quantities. One of these methods 
starts from the distribution differential identity 

V 2 i = -4^5(r) (1) 
r 
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that allows one to write, after specializing for the electron-positron pair fliOfl 

(S(n+)) T = [ 5(n+)* 2 T (K)dR= (2) 

-h I ^4Srt + [Vri ln ^ R )] 2 ^ R 

where R = (ri, T2, r + ) is a point in configuration space, and the trial wave function is 
normalized. Although this integral has a well defined value that can be computed sampling it 
is well known that its variance diverges over the same distribution jlj]. This fact implies that no 
error bound (i.e. standard deviation) can be associated to its value, a dangerous situation one would 
like to avoid. 

Langfelder et al. [ p"4| proposed a possible way to circumvent this problem based on a modified 
importance sampling transformation where V r ?+ i s sampled instead of ^>?r. 

Always starting from Eq. |^, one could also exploit the approach proposed by Assaraf and Caffarel 
]l5[ ] to compute the expectation value needed to obtain nuclear forces by means of the Hellman- 
Feynmann theorem. They showed that a judicious choice of a renormalizcd operator, whose mean 
value is equal to the original one, can reduce the infinite variance to a finite value [ fl6| |. 

A completely different approach was pursued by Alexander and Coldwell |L7|. They proposed 
to compute all the mean values sampling an analytically normalizable distribution function ,g(R), 
so that the normalization Nt of is easily estimated by means of 

M 

w =Af- 1 ^^(R I )/ 5 (R0 (3) 

T i=l 

where the M points sample the normalized g. If a second normalizable distribution g c (R), con- 
strained on the subspace r + = ri, is employed to guide the simulation and to compute Nj, in Eq. 

[| then (8(r |_)) can be easily estimated by the (N^/Nt) 2 ratio. 

Although these three methods represent a step towards the solution of this complicate problem 
in the VMC framework and are currently used for ordinary electronic compounds with success, the 
situation still remains far from being satisfactory for positronic systems. For these systems beyond 
the problem of the method used to compute (8), there is another difficulty: as far as we know, nobody 
has been able to optimize an accurate \&t for a positronic system with more than four electrons. 
More specifically, for large systems explicit correlation between the electrons and the positron has 
been found difficult to introduce This means that the "pile up" of the electron density over the 

positron is not correctly described, therefore giving rise to too small annihilation rates Possible 
sources of this unwanted outcome are the lack of knowledge about the complicate analytical form 
that such an accurate wave function should have, and some drawbacks of the optimization method 
used @. 

In order to go beyond these difficulties, the diffusion Monte Carlo (DMC) method is usually 
employed to sample ^t^q p(J) H)- This technique is able to project out the contribution of 
the excited states from the starting ^>t, allowing the exact calculation of the ground state energy. 

Unfortunately, the 8(r |_) operator does not commute with the Hamiltonian of the system, so that 

the simulation results are only an approximation to the exact mean value when computed by means 
of the mixed estimator 

(S(r-+)) M = J (5(r_ + )* T * rfR (4) 

Although, this value represents a more accurate estimate of the exact (8(r [_)) than (8(r \-))t, it 

has been found that the quality of the results strongly depends on how accurately mimics the 
correct interparticle distributions. 
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Whereas both the spherical well and Gaussian method can be employed to estimate {6(r \-))m 

in Eq. ||, Jiang and Schrader || pointed out that the use of the differential identity Eq. [I] in a 
DMC simulation requires some uncontrolled approximation, since ^t^o is not known analytically 
but only sampled. 

Nevertheless, it has been shown [|l0| that an accurate estimate of (5(r |_)) can be obtained simply 

substituting vtip in Eq. |^ with ^>t^o, if correctly describes the positron-electron distribution. 

A possible solution to the difficulty that DMC meets in estimating the exact expectation values 

is represented by sampling \&q instead of ^t^o, and computing (S(r (_)) without resorting to ^>t 

in any way. This idea rules out the possibility to use Eq. |l|, since one just samples the exact \l/g 
distribution and no analytical information is available about its form. 

In order to overcome this problem, Langfelder et al. [Q proposed to correct (5(r))x by accu- 
mulating the walker weights in a small sphere around the coalescence point. Although this may look 
as a promising way we noticed in our work on positron complexes [^3| that long decaying times 
are needed in order to project out all the excited state contributions and to correctly build the "pile 
up" of the electron density over the positron if poorly describes this feature. This fact produces 
large fluctuations in the weight values, therefore increasing the statistical noise of the results. 

A better approach may be represented by the use of the tagging algorithm proposed by Barnett 
et al. [2J] in connection with the branching step usually employed in DMC. Here, the ratio ^q/^t, 
needed to sample "Jq, is computed by means of the number of daughters of each configuration. 

Moreover, Baroni and Moroni p5[ have recently proposed an alternative algorithm that appears 
to be well suited for this task. This is based on a "Path Integral" view of the DMC algorithm, where 
the branching step has been substituted by an accept/reject step in order to exactly sample 

Unfortunately, these approaches do not solve the problem of the scarce sampling in the volume 
around r = 0, a problem that is present even for small simulation time steps. As stated previously, 
this comes from the finite length and discrete nature of the QMC simulations. As an attempt to 
overcome this difficulty, Langfelder et al. ]l4[ ] implemented in their algorithm the correct sampling of 
the electron-nucleus cusp region as proposed by Umrigar et al. p6| : this, however, does not appear 
straightforward to adopt to correct the sampling of both the electron-electron and electron-positron 
cusps. 

Keeping in mind all the aforementioned problems in estimating (S), we believe the Monte Carlo 
practitioners are left only with the hope of devising an approximate, but hopefully solid and accurate, 
method to compute this observable. 

The main aim of this paper is to discuss and test the accuracy of computing (6) using some 
simple methods based only on the sampling of the positron-electron distribution function without 
any usage of the differential identity Eq. [l|. These methods will be compared with the Gaussian 
approximation discussing relative merits and applicability. Moreover, we apply them to the realistic 
c + LiH model case in order to study the annihilation rate as a function of the internuclear distance R. 
The T2 7 versus R results will allow to compute the vibrationally averaged annihilation rate for this 
system and to discuss molecular environment effects on the annihilation rate itself and on contact 
distribution functions. 

The outlines of this work follow. In Section II we present the basis of the methods. Section III 
describes their applications to model systems for which the exact (6)'s are known. As application 
of this technique, we deal in Section IV with the model e + LiH. Our conclusions and proposals of 
future work are then presented in Section V. 

II. METHODS 

Since we want to develop a method that can be applied to any QMC technique, henceforth we 
will use /(R) to indicate cumulatively *y(R), ^t(R) v I / o(R), or *q(R). Here, R = (ri,r 2 , ...,r + ) is 
a point in configuration space, and r + being respectively the i-th electron and positron positions. 

We are interested in computing the expectation value (S(r |_)) over the distribution /(R), i.e. 
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J f(R)S(r. + )dR 
{S{r - +)) = J /(R)rfR (5) 

Recalling that /(R) is symmetric under any exchange between electrons, Eq. can rewritten as 

{S{r - +)) = f P(r-,r + )dr- d r + < 6 > 

where p(r_,r + ) = N e i e j /(R)dr2...drjv- Introducing the new coordinates R |_ = r + r + and 

r |_ = r + r , after integration over R |_ and spherically averaging over r |_ one gets 



/^(r_ + )^(r_ + )dr_ + / »(r_ + )J(r_ + )r 2 . + dr_ + O(0) 
1 1 /r>(r_+)dr_ + /n(r_+)r?. + dr_ + / n{r^+)r 2 _ + dr_+ 1 j 

where f2(r |_) is the spherically averaged positron-electron distribution. Although these manipu- 
lations are quite straightforward, they highlight that in order to compute (S(r |_)) one must have 

accurate values for both 0(0) and the denominator J £l(r y-)r 2 _ + dr Therefore, both the coales- 
cence region and the tail of the distribution must be correctly described. 

In order to thoroughly present the complexity of the faced problem, Figure 1 shows a typical 
behavior of 0(r |_) as sampled from the model wave function for one electron and one positron 



*i(r_, r+) = exp[-r_ - 0.25r + - 0.25r_+] (8) 

by means of a standard VMC simulation using the Langevin algorithm and the accept /reject step 
Q. This simulation was carried out sampling a grand total of 3.75 x 10 9 configurations and using 
a time step of 0.01 hartree -1 , a fairly small time step for this simple wave function. The sampled 

distribution of r |_ = r was collected on a grid with a bin width of Sr =0.025 bohr, and then the 

number of times r |_ was found inside a given bin was divided by the volume of the spherical crown 

[r — 5r,r], V(r, Sr) = ^-[3r 2 Sr — 3rSr 2 + Sr 3 ]. The so obtained values were attributed to the mean 
radius of the spherical crown, r = 7T(5r[4r 3 — 6r 2 Sr + ArSr 2 — 5r 3 ]/V(r, Sr). This is equivalent to 
approximating il(r) as a straight line inside each spherical crown, a fairly good approximation in 
such a small bin. 

We want to stress that the shape of the distribution was found independent of the time step over 
a broad range of values. This ensures that no systematic bias is present due to the finite time step. 
From Figure 1 one can easily notice the abrupt decrease of the distribution in the region close to 
r = 0. This represents the aforementioned inability of Monte Carlo simulations to correctly sample 
the distribution close to a coalescence point in spite of the large number of sampled configurations. It 
also seems to indicate that, due to this inability, any well-based method (e.g. both the spherical well 
and Gaussian methods) should return an inaccurate answer for smaller that a certain threshold. 
Conversely, all the regions with r > 0.5 bohr seem to be adequately described by the sampled 
distribution, and therefore we propose to analytically continue their shape extrapolating to r = by 
means of a suitable functional form. This idea allows one to exploit the knowledge about the exact 
form of f2 to improve the local description in the small radius regions. For instance, if one samples 
/ = ^q, the exact value of the cusp condition can be used as a way to constrain the model to behave 
correctly. This trick can also be used in both the VMC and DMC simulations, since it is often easy 
to obtain the cusp condition of the sampled / knowing the analytical form of SS*t- This method could 
be implemented in two different ways. First, one could choose an analytical function u(r) to fit £1 
for all the electron-positron distances. This function should be flexible enough to properly describe 
both short range and long range behavior of O. More specifically, close to r = f2(r) behaves like 
exp(— ar), where a is strictly related to the cusp condition. Differently, in the large r regions f2(r) 
follows exp(— (3r), where j3 is dependent on the positron affinity (PA) of the system. A possible 
choice for u>(r) is the Pade-Jastrow form 
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u(r)=N u exp[- ^ ] (9) 
1 + "fr 

where a can be chosen in order to have uj(r) satisfying the correct local behavior close to r = 0. 
The fitted form can successively be used to estimate both f2(0) and the denominator in Eq. |7|. 

Second, if the form of O(r) is more complicated (e.g. it has multiple maxima), it is possible to 
resort to a local fit by Lo(r) in the region close to the cusp in order to get the fl(0) value. Then, the 
normalization integral could be split in two part, one computed using uj(r) and the other directly 
using the sampled distribution. Specifically for the distribution in Figure 1, one could fit the sampled 
VI values in the range [0.5,1.0] bohr constraining oj(r) to have both the exact cusp behavior and to 
have the same value of the sampled Q for r=1.0 bohr. Then, the normalization integral can be 
estimated integrating numerically o>(r) for < r < 1.0, and employing the trapezoidal formula for 
the remaining sampled values. We would like to mention that this necessity is already present for 
small systems like e + Be |f27f . 

Although the two proposed methods are approximate, they might prove themselves to be quite 
accurate in practice, allowing the Monte Carlo practitioner to easily estimate the collision probability 
between two particles. In turn, this will allow to compute the annihilation rate in positronic systems, 
therefore opening the chance to directly compare with the experimental results. 

In order to show that this is exactly the case, in the next section we present the results obtained 

computing (S(r |-)) for some model systems whose exact values are easily obtained using different 

methods. 

Methods similar to ours, although quite different in many details, have been applied by Ortiz 
(2§| and by Mairi Fraser |2{J to the case of a positron embedded in the jellium. Their methods 
were somehow tailored to the specific system under study, so that no direct comparison with our 
proposals can be made. Nevertheless, the results they extracted from the simulations can be useful 
to seek for possible correlations between the magnitude of the "pile up" effect and the local electron 
density. 



III. TEST OF THE METHODS USING MODEL SYSTEMS 

To test the accuracy of the proposed methods, we computed the (S(r h )) expectation value for 

simple model systems composed only by one electron and one positron. More specifically, some 
model wave functions were chosen in order to represent the variety of electron, positron and 
electron-positron distributions that could be found in a positron atomic system. Then, the \&?'s 
were sampled by means of VMC simulations similar to the one discussed above, in order to collect 
the electron-positron distribution Q. 

We selected three model systems as representative of a fairly large class of positron complexes. 
The first one is given by the wave function (see Eq. ||). The second has the analytical wave 
function 

0.15r+- 0.5r2 

* 2 (r_,r+) =exp[-r_ + — + - 0.5r_+] (10) 

where the simple exponential in r + of is substituted by a Pade-Jastrow and the exact cusp 
condition between electron and positron has been introduced. 

To mimic the presence of core and valence shells, we chose as a third function 



r r„+2r 2 n r 15r_-3r 2 nl 0.15r+ - 0.5r 2 
* 3 (r-,r+) = {cxp[ 1 + r _~ ] + 0.001 cxp[ — — ]} exp[ ^ + r+ + - 0.5r_+] (11) 

To compute the exact value (S(r |_)) for these models we used its definition 
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(x(r y _ /^(r-,r + )rf(r_ + )dr + dr_ / tt?(r_,r_)rir- 

+ /*?( r -> r +)*+dr_ r+)dr+dr_ 1 J 

The simple radial integral at the numerator was computed by numerical integration on a grid, while 
the denominator was estimated by means of Eq. | using g(r_ , r+) = A 3 B 3 exp[— 2Ar+ — 2BrJ\/n 2 . 

The two distributions sampled as a function of r from \&2 and ^3 turned out to possess a behavior 
quite similar to Therefore, we avoid to show all of them and refer to Figure 1 as a template for 
such distributions. Due to their smoothness, we fitted them using Eq. || over the range 0.3-10 bohr 
constraining the Pade-Jastrow form to have the exact cusp condition, i.e. -0.5, -1.0, and -1.0, for 
the three respectively. To test the correctness of the chosen fitting range, we modified slightly 
the lower limit without finding statistically meaningful differences. Then, the fitted u> was used to 
estimate both fi(0) and the denominator in Eq. 0by means of numerical integration. The computed 
results, shown as (5(r_+))p a <j e 'i ar e presented in Table [j] together with the "exact" values computed 
using Eq. [l2|. 

During the simulations we also computed the mean values (G(r ^,7)), where G(r (.,7) = 

iV 7 exp[— r?. + /7] is a normalized Gaussian function, for 7 = 0.01, 0.0033, 0.002, and 0.001. This 
technique was proposed by Kenny et al. |30| and successively applied in Ref. and Ref. pT| . 

The (G(r |_, 7)) values were extrapolated to 7 = by fitting them with the simple function a^+b, 

the extrapolation law deduced in Ref. using model systems. The fitting was quite accurate for 
all the three cases, and the results for b = (G(r \-, 0)) are also shown in Table |l|. 

Comparing (S(r — \-))p a de' with the exact results it strikes the really good agreement between 
these two sets of values, the relative error being 1% at most for all the models. It must be pointed 
out that this level of relative accuracy is sufficient to thoroughly compare with the experimental 
data. It is also worth to stress that the application to other model systems gave a similar or better 
relative accuracy, therefore showing the wide applicability of the method. 

As already pointed out previously [|[ |J, also the extrapolated (G(r |_,0)) values are in good 

agreement with the exact results. Nevertheless, one should expect to obtain really inaccurate ap- 
proximations to the exact results using (G(r |_, 7)) with 7 smaller than some threshold value. This 

simple idea is based on the incorrect sampling of the density close to r = shown in Fig. 1, so the 
good agreement found in this and previous works calls for an explanation. This is easily obtained 
superimposing G(r H7) r -+ to the sampled Q(r |_), i.e. comparing the behavior of the two fac- 
tors that form the function whose integrals must be estimated. It turns out that G(r H7) r -+ f° r 

7 > 0.001 has its largest values where ft(r 1-) still behaves correctly, therefore allowing a correct 

estimate of the integrals. Tests carried out using smaller values of 7 gave much worse results than 
the ones reported, so we believe it is safer to limit the values of this parameter to the range 0.001- 
0.01 in order to obtain a meaningful extrapolation. Although this may look some way problematic 

due to the aforementioned difficulties, from the results in Table [l] (G(r (-,0)) appears to be a good 

first estimate of the exact (S(r |_)). In conclusion, we suggest Monte Carlo practitioners to carry 

out always both estimations, i.e. extrapolating (G(r 1^,7)) and fitting the sampled Q, as a way to 

safely estimate (S(r |_)). 

As far as diffusion Monte Carlo and the exact sampling of 'Sq are concerned, the application of 
these methods is straightforward, and no more complications are expected than in the VMC case. 



IV. THE e+LiH SYSTEM 

Having verified the accuracy of the proposed method in computing Dirac's delta mean values, we 
applied it to the calculation of the annihilation rate T2 7 of e + LiH for various internuclear distances 
R. 

Although this system has already been carefully studied employing both QMC methods O, [H| 
and Explicitly Correlated Gaussian (ECG) functions j(| [?], ^] , a description of T 2l as a function of 

the molecular geometry is still lacking. Up to now, there are only (S(r |_)) = (5(r 1+ )) + (<5(r 2 +)) + 

(S(r 3+ )) + (<5(r 4+ )) results at R=3.015 bohr (0.0240(8) from DMC simulations (Sj and 0.024992 
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form ECG calculations ||), at the estimated equilibrium distance R=3.348 bohr (0.027252) 0, 
and the non-adiabatic results of Mitroy and Ryhzikh, 0.034016 and 0.032588 p2[ . These last values 
were obtained using ECG in connection , respectively, with the Stochastic Variational Minimization 
(SVM) and the Frozen-Core SVM (FCSVM) methods, and are roughly 15-20% larger than the ECG 

[F] and DMC [[Tl| clamped-nuclei ones. This unexpected result lead Strasburger || to consider 
the possibility of the flattening of the potential energy curve of e + LiH with respect to the LiH one, a 
feature that may allow the positronic molecule to visit the large internuclear distance region where 
the Dirac's delta mean value is expected to be larger. However, Mitroy and Ryhzikh |?2| pointed 
out that simple ECG's may not represent the best basis functions to describe vibrational nuclear 
motion, and that their vibrational averaged nuclear distances are probably too large. 

In a previous work |3l}| , we computed the complete curve using the DMC technique, showing 
that the flattening is indeed present, and that a strong red shift of the vibrational spectrum with 
respect to LiH must be expected. Unfortunately, due to the high computational cost of our highly 
correlated trial wave functions, at that time we did not compute the behavior of (<5) as a function 
of R. In this work we still adopt the BO approximation, so we predict the annihilation rate for each 
vibrational state of the system and compare our values with the non-adiabatic results of Mitroy and 
Ryhzikh || . 

In order to study the effect of the molecular geometry on the annihilation rate without using the 
computationally expensive wave function used in Ref. [|l) , we decided to employ a model potential 
approach to eliminate the "core" electrons of the Li + fragment. We believe such an approach to 
be physically well grounded, as explained by the following supporting reasons. First, the ECG 
calculations carried out by Strasburger || ^ on the e + LiH system show that the annihilation 
takes place primarily with the two electrons that may be attributed to the H fragment. Secondly, 
the frozen core approximation developed by Mitroy and Ryhzikh Q has been found to describe 
accurately the annihilation process in e + Li, e + Be, LiPs, and e + He 3 S when compared with the 
corresponding all-electron calculations. 

In order to reduce the number of active electrons, we used for e + LiH the model Hamiltonian 



n mod = - \ [V 2 + V 2 + V* ] + V* lod (r i ) + V£ orf (r 2 ) - + — + — L + y + (r+ ) 

2 r H \ r H 2 ri2 r H + r x+ r 2 + 

(13) 

Here, the are interparticle distances, 1 and 2 being the electrons, + the positron, and H the 
hydrogen nucleus. Moreover, the Bardsley's model potential |j| 

_ -l + 10ex p[ -2.202 r , L ,] 

where r^i is the distance between the i-th electron and the Li nucleus, has been used to represent 
the Is 2 Li + core electrons. To model the interaction of the positron with the frozen Li + fragment 
we simply added to the repulsive Coulomb potential of the nucleus the potential of the two frozen 
core electrons as described by an STO-ls orbital with Z=3, obtaining 

0+) = ^ + 6exp[-6r u+ ] -2l^t^±l (15) 

' Lz+ ' Li+ 

where rLi+ is the Li-positron distance. 

To test the accuracy of this model potential, we computed the energy for the ground state of the 
three systems Li, Li~, and LiPs. The energy values are respectively -0.1953(2) hartree, -0.2191(3) 
hartree, and -0.4485(3) hartree. They give an electron affinity of 0.0238(4) hartree, a positron affinity 
(PA) of 0.2294(4) hartree, and a Ps binding energy (BE) of 0.0032(4) hartree. While the electron 
affinity turns out to be in fair agreement with the experimental value, namely 0.023 hartree f3~4j , 
both PA and BE are roughly 0.009 hartree smaller than the best estimate These discrepancies 
might be due to the absence of polarization effects of the core electrons due to the two active 
electrons and the positron, and to a relatively inaccurate representation of the Is 2 electron density 
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by means of a single STO-ls function. Nevertheless, since we are primarily interested in obtaining 
a semiquantitative description of the changes in T2 7 for this system, we believe the approximations 
introduced in the model Hamiltonian to be small enough to allow for a correct prediction of the 
trend for this important observable. 

In order to accurately describe the wave function of the three active particles at the VMC level, 
we employed a trial wave function form similar to the one used in Ref. [ p"l| , but slightly modified 
to include the polarization of the electron and positron density of the PsH fragment due to the 
Li + potential. Specifically, the analytical form in Eq. 7 of Ref. has been multiplied by a 

Pade-Jastrow factor depending on the z coordinate of each particle: the z axis was chosen as the 
LiH bond axis, the H nucleus being located at the origin and the Li on the negative z axis. The 
wave function parameters were fully optimized for every nuclear distance R minimizing the variance 
of the local energy over a fixed sample of configurations |3^, This procedure is already well 
described in the literature Jl3|], so we skip the unnecessary details. The ensemble of walkers used 
in the optimization was generated by DMC simulations in order to bias the walkers distribution 
towards the exact density. Usually, four or five optimization steps were carried out for each R. 

We started the optimization process of the wave function at R=20 plf . At the end of the 
optimization procedure, R was decreased and the wave function re-optimized for the new distance. 
This procedure gives the chance to monitor the changes of the wave function with R, but might 
increase the possibility to remain stuck in a local minimum in the parameter space during the 
optimization. 

Having optimized at VMC level the approximate wave functions for various distances, these were 
employed in long DMC simulations to project out the remaining excited state contributions and 
to compute more accurate mixed expectation values. For all the simulations, a time step of 0.005 
hartree -1 was used, together with a population of 9000 walkers. These two simulation parameters 
were found adequate to make statistically negligible both the time step bias and the population 
effect in the DMC simulations. 

The DMC results for the energy and for the (6(r |_)) of this model system are shown in Ta- 
ble |. There, the energy values represent the ground state energy of the model Hamiltonian Eq. 
fl3| , (<5(r \-))p a de' ar e the total collision probabilities estimated using the electron-positron distri- 
bution, i.e. (5(ri+)) + (S(r2+)), while (G(r |-,0)) are the extrapolated Gaussian values. Here, the 

electron-positron distributions were fitted with the function in Eq. |p, constraining its cusp to be 
-0.5+cusp(^ T ). 

The energy values obtained in Ref. |3l), after having subtracted the repulsion 1/R between 
the H nucleus and the Li + core and the total energy of the Li + fragment (-7.279913 hartree ||) 
to estimate the leptonic energy of the PsH moiety, are shown in Figure 2 together with the DMC 
results obtained in this work. 

For R>4, the results from the model system follow closely the more accurate all-electron FN- 
DMC values, showing that the model potential correctly describes the polarization of PsH due to the 
interaction with Li + . For shorter distances the approximation of considering Li + frozen is no longer 
accurate, so that a discrepancy between the two sets of results is expected. It is also important to 
remember that the Bardsley's potential was tailored only to describe the atom in its ground and 
valence excited states, not to describe correctly bonds in molecules. 

Figure 3 shows the two computed (5(r |_)) values as a function of the intcrnuclcar distance R, 

together with the ECG results of Strasburger §, § and the DMC result of Ref. |l[. These last 
results can be used to evaluate the total accuracy of our computed collision probabilities. 

From Fig. 3 it is clear that both the (S(r \.))p a de' and the (G(r ^,0)) are in good agreement, 

the second differing from the first one by two standard deviations at most. This allows us to think 
that we are accurately estimating the mixed distribution mean value obtained by the standard DMC 
technique. Improved results could be obtained only by sampling \&q. 

As far as the total accuracy is concerned, at R=3.015 bohr both all-electron ECG Q and 
DMC results jll| appear to be smaller than the model ones by roughly 7%. Instead, on going 
towards large R the mean value seems to correctly converge to the very accurate ECG value, namely 
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(5(r |-)) p S H =0.04874 ||, and to the new DMC 0.0486 estimate carried out in this work using the 

electron-positron distribution. From this comparison, one could expect our estimation to slightly 
degrade on going towards small R, without becoming embarrassing inaccurate to create concerns 
about the usefulness of this model system. 

The overall trend of the collision probability shows a net decrease on going towards short R, a 
feature already suggested by Mitroy and Ryzhikh [ [32| . This can be easily understood remembering 
that the Li + model potential repels the positron, while attracting the two electrons. It is not easy 
to infer any possible analytical model to describe these joined effects, although for large R one could 
propose a limiting 1 /R 2 form due to the polarization of the two distributions by the electric field of 
Li + . We show in the Appendix that this reasoning is indeed correct for any observable by means 
of first order Perturbation Theory. An "experimental" evidence that this is the case is given by the 
fairly good fitting of the (5) results at R=10, 15, and 20 bohr with the simple form 0.0486 + b/R 2 , 
where b = -0.20448. 

Various other mean values were computed during the DMC simulations in order to obtain some 
physical insight on the electron and positron density behavior. Figure 4 reports the mean value 
of the z coordinate for the two particles giving information on the polarization of the two lepton 
densities. It is clear that the positronic distribution is polarized by the model potential in the 
opposite direction to the electronic one. Moreover, it appears to be more easily polarized than the 
electronic one always showing larger (z) values. This fact can be easily explained by noticing that 
the positron distribution is more diffuse that the electron one, so that it is more strongly repelled by 
the electric field. Interestingly, at R=2 and 2.5 bohr the electron distribution reverse its polarization 
showing (z) > 0. We believe this effect to be due to the repulsive region of the Li + core potential 
that pushes away the electrons for such small nuclear distances, therefore modeling the exchange 
effect created by an antisymmetric wave function. 

Figure 5 shows the average values for the electron-electron and electron-positron distances. While 
the electron-positron mean distance (r |_) increases monotonically upon decreasing R, the electron- 
electron distance (r ) shows a shallow maximum around R=8 bohr and a deep minimum around 

R=2.5 bohr. We believe the maximum to be due to the competition between the positron and the 
Li + model potential to bind an electron. More specifically, although polarized towards positive z, 
the positron still attracts one of the two electrons to form the Ps sub-cluster, increasing the distance 
from the second electron that is free to be polarized in the direction of the Li + core. On going 
towards smaller R, the positron is pushed far out the bond region, loosing its ability to polarize the 
electrons that are now both strongly attracted by the model potential. This interaction leads them 
to move in the small volume between H and Li + , therefore decreasing their mean distance. Then, for 
R smaller than 2.5 bohr, the electron-core repulsion pushes the electrons outwards the bond region, 
with the net effect to increase their mean distance. This effect has also been observed by plotting 
the intracule electron distributions obtained during the DMC simulations. 

It is worth to mention that similar conclusions can be drawn analysing the VMC results obtained 
as a by-product of the optimization stages. 

Having studied the overall behavior of (S) , (z) , and (r) , we now turn to compute the vibrationally 
averaged annihilation probabilities. To obtain these quantities, we interpolated our (S) results by 
means of the analytical form A(R) = 0.0486 - 2aR/{\ + bR + cR 2 + dR 3 ). The fitted parameters 
are a = 1.05945, b = 97.5779, c = -37.9705, and d = 12.2715. Then, the potential energy curve of 
e + LiH obtained in Ref. ]3lJ was fitted with the modified Morse potential 

V M {R) = -8.0699 + A{1 - cxp[-B(R - C*)]} 2 - A — D{1 — exp[-(i?/ J F)] 6 }/(2i? 4 ) (16) 

obtaining A = 0.03444 hartree, B = 0.72030 bohr" 1 , C = 3.3060 bohr, D = 21.1796 bohr" 3 , and 
F = 5.88217 bohr. The last term in Eq. |l6| has been introduced in order to correctly represent 
the charge-induced dipole interaction between Li + and PsH. The nuclear Schrodinger equation for 
this potential was then solved using the grid method proposed by Tobin and Hinze ]37| , and the 
numerical wave functions <f> v (R) were then used to compute vibrationally averaged mean values for 
zero total angular momentum. More specifically, we computed 
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_ fdR4l(R)Q(R) 
{U} »- fdRft(R) (17) 

where O(R) is A(R) or any other function of R. 

In Table || we show the results for (S) u and (R) v computed over the first 16 bound vibrational 
states. The (£)„ values increase in an almost linear fashion on going towards large v, as expected 
by the steady increase of (R) v due to the vibrational excitation. Comparing (8)q =0.0295 with the 
value of A(i?) at the equilibrium distance of our fitted potential, namely 0.0291 at 3.353 bohr, it 
appears that the ground level vibrational motion only slightly increases the probability of collision 
between the electrons and the positron with respect to the one at the equilibrium distance. This 
finding is in line with the small difference between the equilibrium distance and the average nuclear 
distance Rq =3.42 bohr. We relate these outcomes to the almost linear behavior of (A(R)), and to 
the shape of i? 2 0§(i?) in the region around the potential minimum, where it resembles a Gaussian. 

Although our vibrationally averaged result for v = (S)o= 0.0295 appears to be roughly 8% 
larger than the ECG result (0.027252) at the equilibrium distance 3.348 bohr suggesting a 
fairly large effect of the nuclear motion, we believe this outcome is primarily due to the 7% larger 
collision probabilities computed using our model system. These two evidences seem to rule out the 
Strasburger's suggestion || of a large increase of the collision probability due to the quantum nuclear 
motion for the v = state. They indicate that approximating the averaged collision probability for 
the vibrational ground state simply by using its value at the equilibrium distance could be a fairly 
accurate procedure. Moreover, these conclusions agree with Mitroy and Ryhzikh p2| warnings 
that both the SVM and FCSVM results, although proving the overall stability of e + LiH, are not 
well converged to the exact ones. For instance, their (R)o values, respectively 4.182 bohr and 3.964 
bohr, are larger than the minimum of the ECG and DMC potential curves by more than 0.5 bohr. 
This discrepancy cannot be accounted for by the zero point motion of the positron complex. These 
larger distances between the two fragments Li + and PsH in the non-adiabatic treatment imply a 
reduced distortion of the lepton densities of PsH with respect to the one at the Born-Oppcnheimer 
equilibrium, and therefore a too large annihilation rate. However, it is interesting to notice that 
both FCSVM (0.032588) and SVM (0.034016) §§ collision probabilities are really close to our 
Born-Oppenhcimer one at R =4.0 bohr. In our view, this agreement stresses, again, the importance 
of the local electric field in defining the collision probability and the overall accuracy of the SVM 
approach in describing the relative densities an a positronic complex. 

As far as the behavior of (6(R)) V is concerned, the steady increase on going towards large v 
indicates that the annihilation rate does depend on the quantum vibrational state of the molecule. 
Although the trend of these results could be specific of the e + LiH system and perhaps of other 
polar molecules as well, it strongly indicates that any theory formulated to describe "on the fly" 
annihilation of e + due to Feshbach resonances must include this effect in order to go beyond "order 
of magnitude comparison" fl38|| and to predict accurately the annihilation rate. In our view, this 
opens a new avenue of exploration in positron physical chemistry where the understanding of the 
vibrational motion effect on positron annihilation by molecular systems is of prime importance. 



V. CONCLUSIONS 

In this work we have critically compared methods that may be useful to compute the annihilation 
rate in positronic systems in the framework of the QMC methods. Moreover, we have presented 
a simple, but nevertheless solid and accurate, method based only on the interparticle distribution 
sampling. After having tested it using model systems, we employed the method to compute (<5(r_ + )) 
for e + LiH for several internuclear distances. These results allowed us to discuss many interesting 
features of this positronic complex, and to predict that the annihilation probability increases upon 

increasing the vibrational quantum number v. We notice that a similar behavior of (6(r |_)) may 

be expected also for e + LiF due to the polarization of the positronic density of the PsF fragment by 
the Li + core. The situation could be quite different for the e + BeO case where the positron density is 
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expected to be centered on Be at large nuclear distances (the two fragments e + Be and O have lower 
total energy than Be + and PsO ]39|]), and to move on the O side of the molecule when the distance 
decreases. This effect is due to the electron transfer from Be to O that creates the large molecular 

dipole moment. From our experience on these systems we expect (S(r |_)) for e + Be to be smaller 

than the one for the polar molecule, so the vibrationally excited states close to the dissociation 
threshold may have smaller annihilation rates than the ground vibrational state. 

It is also interesting to speculate on the behavior of the annihilation rate versus the vibrational 
quantum number for other simple systems like e + Li2 and e + Be2- Here, the symmetry of the systems 
can play an important role in defining the annihilation rate. For instance, decreasing the nuclear 
distance one may expect to find the positron localized between the two atomic fragments due to its 
ability to polarize the two atomic electron densities: in this situation the annihilation rate could be 
quite different from the atomic one. Although it is easy to infer the existence of a bound state for 
these complexes at large nuclear distances employing the basic Valence Bond resonance idea 

e+A 1 +A 2 < >A 1 +e+A 2 (18) 

it still remains to demonstrate the stability of these systems for nuclear distances close to the 
equilibrium geometry of the neutral parent molecules. 

As rule of thumb to predict the stability of a non-polar molecule, one can use the adiabatic 
ionization potential (AFP) as proposed by Mitroy et al. f4Cj| . For the X ground state of Li2 the 
experimental AIP is 0.189 hartree Q], slightly lower than the atomic one, 0.19814 hartree E |. 
Also for the X 1 S+ ground state of Be2, one might expect a similar lowering of the AIP with respect 
to the atomic one, 0.343 hartree Q, so that a value of around 0.335 hartree could be regarded as 
a safe upper bound to the true AIP. Both these values fall inside the upper and lower IP limits for 
positron binding obtained by Mitroy et al. |Q for one and two valence electron atoms, therefore 
suggesting that the two complexes should be stable. We understand that this model is just a rough 
approximation for our molecular systems p3") . Nevertheless, a positron bound to an atom or a 
molecule is always characterized by a quite diffuse density. This allows one to neglect some of the 
real features of the electron density close to the nuclei as a first approximation, and focus only on the 
asymptotic properties of the positron cloud that are correlated to the IP and to the polarizability. 

As far as Be 2 is concerned, the AIP larger than the Ps binding energy (0.25 hartree) suggests 
a mechanism based on the electron cloud polarization as responsible of the binding. Moreover, the 
lowering of the AIP with respect to the atomic one, and the large polarizability of this molecule 
(roughly twice the atomic one) seem to indicate its ability to form a stronger bond with the positron 
than the Be atom alone [Q. Conversely, Li 2 has an AIP smaller than the Ps binding energy, 
suggesting that the polarization of the Ps cluster may be held responsible for the positron binding. 
However, Li is close to the lower stability threshold of the positron-atom complexes, and we do 
not feel confident in proposing the stability of the molecular complex with respect to the Ps+LiJ 
dissociation pathway. 

Right now, QMC methods are the best suited computational techniques to carry out such a 
study since 6 and 8 electron systems are too large to be studied with ECG's unless the frozen core 
approximation is used. With the addition to the QMC "bag of tricks" of a robust method for 
computing annihilation rates such a study could become routine in molecular physical chemistry, 
allowing the exploration of many interesting features of these "exotic" compounds. 

Moreover, many more other technically oriented applications could be devised. Positronium 
annihilation in polymers and membranes, positron annihilation in silicon nanocluster, nanodevices, 
fullerenes and carbon nanotubes are just few of them that could be quite easily interpreted with the 
help of such a method. 

Our hope is that this work will help this kind of applications to blossom and to lead to a better 
understanding of the basic interaction schemes that positron has with ordinary matter. 
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APPENDIX 



In this Appendix we show that the correction to the expectation value of every observable O for 
PsH interacting with the Li + core follows the limiting analytical form R~ 2 for large R. This is indeed 
a general results for positronic atom systems immersed in a weak electric field. From Perturbation 
Theory one can write the first order corrected wave function for the ground state as 

*P = *o 0) + E £ S^r *™ = ^ + £ ^ (19) 

where \l/^ ' ) are the eigenstates of the unperturbed Hamiltonian (i.e. the PsH), E^p its eigenvalues, 
and V the perturbation potential. The expectation value of the observable O can be computed using 

c*Oc 

where Oji = J O^^dR., while c = {1, cf \ ■■■}■ If the perturbation potential is small with 

respect to the total energy, then c^ "* <§; 1 and c'lc = 1 + J2i=£o( c i°' ') 2 — 1- This fact allows to 
introduce a further approximation into Eq. ^0], 

(O) ~ c*Oc=<0>oo + 2]T Cl (0) (O) 0s :+ ]T cf cf ) (0) Jl ^(0) o + 2^cf (0)o, (21) 

showing that the first order change in the expectation value (O) is linearly dependent on the 's. 

For our specific case, namely PsH interacting with the Coulomb potential of Li + at distance R, 
for R — > oo the perturbing interaction potential can be written as 

where the molecular geometry is as in the main text, while qu are the leptonic charges. This 
approximation is equal to consider the electron and positron densities constant in a plane parallel 
to the xy plane. Introducing this approximation in the integrals in Eq. 09 one gets 



/ ^ffl^fti / <^M 0) <« + 5> / *l°>%*V>dR= (23) 
J k J k 

R 2 

This result shows that the c, in Eq. |2(] and ^l] are proportional to l/R 2 , therefore proving that 
this is also the analytical form of the leading correction to the unperturbed ground state expectation 
values O o- 
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Figure Captions: 



Figure 1: Behavior of f2(r (_) sampled from $1. 

Figure 2: Energy of the PsH moiety computed from Eq. 13 and from the e + LiH results of Ref. 
PH l as explained in the text. 

Figure 3: Computed (S) results for c + LiH system at various internuclear distances. 

Figure 4: Electron and positron mean value of the z coordinate as a function of the internuclear 
distance. 

Figure 5: Electron-electron and electron-positron mean distances as function of the internuclear 
distance. 
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( r -+)) exact 


{5{r- + ))pade' 


7 


(G(r_+,7)) 




0.022602(7) 


0.02229 


0.0100 
0.0033 
0.0020 
0.0010 
0.0000 


0.02140(8) 
0.0219(1) 
0.0222(2) 
0.0223(3) 
0.0228(3) 




0.10981(1) 


0.11052 


0.0100 
0.0033 
0.0020 
0.0010 
0.0000 


0.0979(2) 
0.1028(3) 
0.1043(4) 
0.1058(7) 
0.1095(7) 








0.0100 
0.0033 
0.0020 


0.0842(2) 
0.0886(3) 
0.0901(4) 




0.09399(1) 


0.09382 


0.0010 
0.0000 


0.0916(7) 
0.0950(7) 



Table 1: (S(r |_)) expectation values for the three model systems VPi, \&2, and ^3. The "exact" 

values are computed by Eq. [l^. (S(r |_)) p a de' are computed fitting Eq. [|to the sampled distribution 

as explained in the text. 7 is the width of the Gaussian used to compute (G(r |_, 7)). 
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R 


(E) 




<G(r_+,0)> 


2.0 


-1.16437(9) 


0.01884 


0.0182(2) 


2.5 


-1.15936(7) 


0.02178 


0.0210(2) 


3.0 


-1.13196(6) 


0.02684 


0.0256(2) 


3.5 


-1.09798(5) 


0.03000 


0.0288(2) 


4.0 


-1.06485(5) 


0.03334 


0.0324(2) 


6.0 


-0.96985(4) 


0.04166 


0.0402(2) 


8.0 


-0.91995(3) 


0.04544 


0.0436(3) 


10.0 


-0.89160(3) 


0.04666 


0.0448(3) 


15.0 


-0.85634(2) 


0.04754 


0.0458(3) 


20.0 


-0.83932(2) 


0.04794 


0.0468(3) 


oo 


-0.78918(1) 


0.04860 


0.0484(3) 



Table 2: Lepton energies, (S(r |_)) p a de' , and (G(r h , 0)) mean values for the c+LiH model system. 

All quantities in atomic units. 
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.0295 


3 


.423 


1 





,0305 


3 


.571 


2 





.0315 


3 


.732 


3 





.0325 


3 


.906 


4 





.0334 


4 


.094 


5 





.0344 


4 


.294 


6 





.0352 


4 


.503 


7 





.0361 


4 


.718 


8 





.0369 


4 


.938 


9 





.0376 


5 


.166 


10 





.0383 


5 


.408 


11 





.0390 


5 


.675 


12 





.0397 


5 


.985 


13 





.0405 


6 


.356 


14 





.0414 


6 


.803 


15 





.0421 


7 


.243 



Table 3: Vibrationally averaged (S(r y.)) v and (i?)„ values for the e+LiH model system. All quan- 
tities in atomic units. 
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